/*


!
!  Dalton, a molecular electronic structure program
!  Copyright (C) The Dalton Authors (see AUTHORS file for details).
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License version 2.1 as published by the Free Software Foundation.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  Lesser General Public License for more details.
!
!  If a copy of the GNU LGPL v2.1 was not distributed with this
!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
!

!

*/
/*-*-mode: C; c-indentation-style: "bsd"; c-basic-offset: 4; -*-*/
/* fun-revpbe.c:

   Automatically generated code implementing revPBEx functional and
   its derivatives. It is generated by func-codegen.pl being a part of
   a "Automatic code generation framework for analytical functional
   derivative evaluation", Pawel Salek, 2005

    This functional is connected by making following changes:
    1. add "extern Functional revpbeFunctional;" to 'functionals.h'
    2. add "&revpbeFunctional," to 'functionals.c'
    3. add "fun-revpbe.c" to 'Makefile.am', 'Makefile.in' or 'Makefile'.

    This functional has been generated from following input:
    ------ cut here -------
xa:sqrt(grada*grada)/rhoa^(4/3);
xb:sqrt(gradb*gradb)/rhob^(4/3);
 
R:1.245;
d:0.066725;
mu:d*%PI^2/3;
Sa:xa/(2*(6*%PI^2)^(1/3));
Sb:xb/(2*(6*%PI^2)^(1/3));
 
F(S):=1+R-R/(1+mu*S^2/R);
Ea(n):=-3/(4*%PI)*(3*%PI^2)^(1/3)*n^(4/3)*F(Sa);
Eb(n):=-3/(4*%PI)*(3*%PI^2)^(1/3)*n^(4/3)*F(Sb);
 
K(rhoa,grada,rhob,gradb,gradab):=0.5*(Ea(2*rhoa)+Eb(2*rhob));


    ------ cut here -------
*/

 
/* strictly conform to XOPEN ANSI C standard */
#if !defined(SYS_DEC)
/* XOPEN compliance is missing on old Tru64 4.0E Alphas and pow() prototype
 * is not specified. */
#define _XOPEN_SOURCE          500
#define _XOPEN_SOURCE_EXTENDED 1
#endif
#include <math.h>
#include <stddef.h>
#include "general.h"

#define __CVERSION__
 
#include "functionals.h"
 
/* INTERFACE PART */
static integer revpbe_isgga(void) { return 1; } /* FIXME: detect! */
static integer revpbe_read(const char *conf_line);
static real revpbe_energy(const FunDensProp* dp);
static void revpbe_first(FunFirstFuncDrv *ds,   real factor,
                         const FunDensProp* dp);
static void revpbe_second(FunSecondFuncDrv *ds, real factor,
                          const FunDensProp* dp);
static void revpbe_third(FunThirdFuncDrv *ds,   real factor,
                         const FunDensProp* dp);
static void revpbe_fourth(FunFourthFuncDrv *ds,   real factor,
                          const FunDensProp* dp);
 
Functional revPBExFunctional = {
  "revPBEx",       /* name */
  revpbe_isgga,   /* gga-corrected */
   1,
  revpbe_read,
  NULL,
  revpbe_energy,
  revpbe_first,
  revpbe_second,
  revpbe_third,
  revpbe_fourth
};
 
/* IMPLEMENTATION PART */
static integer
revpbe_read(const char *conf_line)
{
    fun_set_hf_weight(0);
    return 1;
}

static real
revpbe_energy(const FunDensProp *dp)
{
    real res;
    real rhoa = dp->rhoa, rhob = dp->rhob;
    real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;

    real t1, t2, t3, t4, t5;

    t1 = pow(2.0,0.333333333333333);
    t2 = pow(3.0,0.333333333333333);
    t3 = 1/pow(3.141592653589793,0.333333333333333);
    t4 = 1/pow(6.0,0.666666666666667);
    t5 = pow(3.141592653589793,0.666666666666667);

   /* code */
    res = 0.5*(-1.5*t1*t2*t3*pow(rhob,1.333333333333333)*
        (2.245-1.245/(0.004466198125837*t4*t5*pow(gradb,2.0)/pow(rhob,
        2.666666666666667)+1.0))-1.5*t1*t2*t3*pow(rhoa,1.333333333333333)*
        (2.245-1.245/(0.004466198125837*t4*t5*pow(grada,2.0)/pow(rhoa,
        2.666666666666667)+1.0)));

    return res;
}

static void
revpbe_first_helper(real rhoa, real grada, real *res)
{    real t1, t2, t3, t4, t5, t6, t7;

    t1 = pow(2.0,0.333333333333333);
    t2 = pow(3.0,0.333333333333333);
    t3 = 1/pow(6.0,0.666666666666667);
    t4 = pow(3.141592653589793,0.333333333333333);
    t5 = pow(grada,2.0);
    t6 = 0.004466198125837*pow(3.141592653589793,0.666666666666667)*
        t3*t5/pow(rhoa,2.666666666666667)+1.0;
    t7 = 1/pow(t6,2.0);

   /* code */
    res[0] = 0.5*(0.022241666666667*t1*t2*t3*t4*t5*t7/pow(rhoa,
        2.333333333333333)-2.0*t1*t2*(2.245-1.245/t6)*pow(rhoa,0.333333333333333)/
        t4);
    res[1] = -0.008340625*t1*t2*t3*t4*t7*grada/pow(rhoa,1.333333333333333);
}

static void
revpbe_first(FunFirstFuncDrv *ds, real factor, const FunDensProp *dp)
{
    real res[2];

    revpbe_first_helper(dp->rhoa, dp->grada, res);
   /* Final assignment */
    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        revpbe_first_helper(dp->rhob, dp->gradb, res);
    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

}

static void
revpbe_second_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13;

    t1 = pow(2.0,0.333333333333333);
    t2 = pow(3.0,0.333333333333333);
    t3 = 1/pow(6.0,0.666666666666667);
    t4 = pow(3.141592653589793,0.333333333333333);
    t5 = pow(grada,2.0);
    t6 = 0.004466198125837*pow(3.141592653589793,0.666666666666667)*
        t3*t5/pow(rhoa,2.666666666666667)+1.0;
    t7 = 1/pow(t6,2.0);
    t8 = 1/pow(rhoa,2.333333333333333);
    t9 = 1/t4;
    t10 = 2.245-1.245/t6;
    t11 = 1/pow(rhoa,1.333333333333333);
    t12 = 1/pow(6.0,0.333333333333333);
    t13 = 1/pow(t6,3.0);

   /* code */
    res[0] = 0.5*(0.022241666666667*t1*t2*t3*t4*t5*t7*t8-
        2.0*t1*t10*t2*t9*pow(rhoa,0.333333333333333));
    res[1] = -0.008340625*t1*t2*t3*t4*grada*t7*t11;
    res[2] = 0.5*(2.773975767886205E-4*t1*t12*t13*t2*pow(grada,
        4.0)/pow(rhoa,6.0)-0.022241666666667*t1*t2*t3*t4*t5*t7/pow(rhoa,
        3.333333333333333)-0.666666666666667*t1*t10*t2*t9/pow(rhoa,
        0.666666666666667));
    res[3] = 0.5*(0.022241666666667*t1*t2*t3*t4*grada*t7*
        t8-2.0804818259146537E-4*t1*t12*t13*t2*pow(grada,3.0)/pow(rhoa,
        5.0));
    res[4] = 7.801806847179953E-5*t1*t12*t13*t2*t5/pow(rhoa,
        4.0)-0.008340625*t1*t2*t3*t4*t7*t11;

}

static void
revpbe_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[5];
 
    revpbe_second_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        revpbe_second_helper(dp->rhob, dp->gradb, res);
    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

}

static void
revpbe_third_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21;

    t1 = pow(2.0,0.333333333333333);
    t2 = pow(3.0,0.333333333333333);
    t3 = 1/pow(6.0,0.666666666666667);
    t4 = pow(3.141592653589793,0.333333333333333);
    t5 = pow(grada,2.0);
    t6 = 0.004466198125837*pow(3.141592653589793,0.666666666666667)*
        t3*t5/pow(rhoa,2.666666666666667)+1.0;
    t7 = 1/pow(t6,2.0);
    t8 = 1/pow(rhoa,2.333333333333333);
    t9 = 1/t4;
    t10 = 2.245-1.245/t6;
    t11 = 1/pow(rhoa,1.333333333333333);
    t12 = 1/pow(6.0,0.333333333333333);
    t13 = pow(grada,4.0);
    t14 = 1/pow(t6,3.0);
    t15 = 1/pow(rhoa,6.0);
    t16 = 1/pow(rhoa,3.333333333333333);
    t17 = pow(grada,3.0);
    t18 = 1/pow(rhoa,5.0);
    t19 = 1/pow(rhoa,4.0);
    t20 = pow(3.141592653589793,1.666666666666667);
    t21 = 1/pow(t6,4.0);

   /* code */
    res[0] = 0.5*(0.022241666666667*t1*t2*t3*t4*t5*t7*t8-
        2.0*t1*t10*t2*t9*pow(rhoa,0.333333333333333));
    res[1] = -0.008340625*t1*t2*t3*t4*grada*t7*t11;
    res[2] = 0.5*(-0.666666666666667*t1*t10*t2*t9/pow(rhoa,
        0.666666666666667)+2.773975767886205E-4*t1*t12*t13*t14*t15*
        t2-0.022241666666667*t1*t2*t3*t4*t5*t7*t16);
    res[3] = 0.5*(0.022241666666667*t1*t2*t3*t4*grada*t7*
        t8-2.0804818259146537E-4*t1*t12*t14*t17*t18*t2);
    res[4] = 7.801806847179953E-5*t1*t12*t14*t19*t2*t5-0.008340625*
        t1*t2*t3*t4*t7*t11;
    res[5] = 0.5*(5.258108117653037E-7*t1*t2*t20*t21*pow(grada,
        6.0)/pow(rhoa,9.666666666666666)-0.00194178303752*t1*t12*t13*
        t14*t2/pow(rhoa,7.0)+0.084024074074074*t1*t2*t3*t4*t5*t7/pow(rhoa,
        4.333333333333333)+0.444444444444444*t1*t10*t2*t9/pow(rhoa,
        1.666666666666667));
    res[6] = 0.5*(-3.943581088239779E-7*t1*t2*t20*t21*pow(grada,
        5.0)/pow(rhoa,8.666666666666666)+0.001317638489746*t1*t12*
        t14*t15*t17*t2-0.051897222222222*t1*t2*t3*t4*grada*t7*t16);
    res[7] = 0.5*(2.957685816179834E-7*t1*t13*t2*t20*t21/
        pow(rhoa,7.666666666666667)+0.022241666666667*t1*t2*t3*t4*
        t7*t8-8.321927303658615E-4*t1*t12*t14*t18*t2*t5);
    res[8] = 2.340542054153986E-4*t1*t12*t14*t19*t2*grada-
        1.109132181067438E-7*t1*t17*t2*t20*t21/pow(rhoa,6.666666666666667);

}

static void
revpbe_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[9];
 
    revpbe_third_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];

    ds->df3000 += factor*res[5];
    ds->df2010 += factor*res[6];
    ds->df1020 += factor*res[7];
    ds->df0030 += factor*res[8];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        revpbe_third_helper(dp->rhob, dp->gradb, res);

    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

    ds->df0300 += factor*res[5];
    ds->df0201 += factor*res[6];
    ds->df0102 += factor*res[7];
    ds->df0003 += factor*res[8];

}

static void
revpbe_fourth_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30, t31, t32;

    t1 = pow(2.0,0.333333333333333);
    t2 = pow(3.0,0.333333333333333);
    t3 = 1/pow(6.0,0.666666666666667);
    t4 = pow(3.141592653589793,0.333333333333333);
    t5 = pow(grada,2.0);
    t6 = 1/pow(rhoa,2.666666666666667);
    t7 = 0.004466198125837*pow(3.141592653589793,0.666666666666667)*
        t3*t5*t6+1.0;
    t8 = 1/pow(t7,2.0);
    t9 = 1/pow(rhoa,2.333333333333333);
    t10 = 1/t4;
    t11 = 2.245-1.245/t7;
    t12 = 1/pow(rhoa,1.333333333333333);
    t13 = 1/pow(6.0,0.333333333333333);
    t14 = pow(grada,4.0);
    t15 = 1/pow(t7,3.0);
    t16 = 1/pow(rhoa,6.0);
    t17 = 1/pow(rhoa,3.333333333333333);
    t18 = pow(grada,3.0);
    t19 = 1/pow(rhoa,5.0);
    t20 = 1/pow(rhoa,4.0);
    t21 = pow(3.141592653589793,1.666666666666667);
    t22 = pow(grada,6.0);
    t23 = 1/pow(t7,4.0);
    t24 = 1/pow(rhoa,9.666666666666666);
    t25 = 1/pow(rhoa,7.0);
    t26 = 1/pow(rhoa,4.333333333333333);
    t27 = pow(grada,5.0);
    t28 = 1/pow(rhoa,8.666666666666666);
    t29 = 1/pow(rhoa,7.666666666666667);
    t30 = 1/pow(rhoa,6.666666666666667);
    t31 = pow(3.141592653589793,2.333333333333333);
    t32 = 1/pow(t7,5.0);

   /* code */
    res[0] = 0.5*(0.022241666666667*t1*t2*t3*t4*t5*t8*t9-
        2.0*t1*t10*t11*t2*pow(rhoa,0.333333333333333));
    res[1] = -0.008340625*t1*t2*t3*t4*grada*t8*t12;
    res[2] = 0.5*(-0.666666666666667*t1*t10*t11*t2/pow(rhoa,
        0.666666666666667)+2.773975767886205E-4*t1*t13*t14*t15*t16*
        t2-0.022241666666667*t1*t2*t3*t4*t5*t8*t17);
    res[3] = 0.5*(0.022241666666667*t1*t2*t3*t4*grada*t8*
        t9-2.0804818259146537E-4*t1*t13*t15*t18*t19*t2);
    res[4] = 7.801806847179953E-5*t1*t13*t15*t2*t20*t5-0.008340625*
        t1*t2*t3*t4*t8*t12;
    res[5] = 0.5*(0.444444444444444*t1*t10*t11*t2/pow(rhoa,
        1.666666666666667)+0.084024074074074*t1*t2*t3*t4*t5*t8*t26-
        0.00194178303752*t1*t13*t14*t15*t2*t25+5.258108117653037E-7*
        t1*t2*t21*t22*t23*t24);
    res[6] = 0.5*(-3.943581088239779E-7*t1*t2*t21*t27*t23*
        t28+0.001317638489746*t1*t13*t15*t16*t18*t2-0.051897222222222*
        t1*t2*t3*t4*grada*t8*t17);
    res[7] = 0.5*(0.022241666666667*t1*t2*t3*t4*t8*t9-8.321927303658615E-4*
        t1*t13*t15*t19*t2*t5+2.957685816179834E-7*t1*t2*t21*t14*t23*
        t29);
    res[8] = 2.340542054153986E-4*t1*t13*t15*t2*t20*grada-
        1.109132181067438E-7*t1*t2*t21*t18*t23*t30;
    res[9] = 0.5*(2.504933612854253E-8*t1*t2*t3*t31*t32*pow(grada,
        8.0)/pow(rhoa,13.33333333333333)-8.763513529421729E-6*t1*t2*
        t21*t22*t23/pow(rhoa,10.66666666666667)+0.014640427663844*
        t1*t13*t14*t15*t2/pow(rhoa,8.0)-0.370694444444444*t1*t2*t3*
        t4*t5*t8/pow(rhoa,5.333333333333333)-0.740740740740741*t1*
        t10*t11*t2*t6);
    res[10] = 0.5*(-1.87870020964069E-8*t1*t2*t3*t31*t32*
        pow(grada,7.0)/pow(rhoa,12.33333333333333)+0.172990740740741*
        t1*t2*t3*t4*grada*t8*t26-0.008553091950982*t1*t13*t15*t18*
        t2*t25+5.915371632359669E-6*t1*t2*t21*t27*t23*t24);
    res[11] = 0.5*(1.4090251572305177E-8*t1*t2*t22*t3*t31*
        t32/pow(rhoa,11.33333333333333)+0.004438361228618*t1*t13*t15*
        t16*t2*t5-3.844991561033785E-6*t1*t2*t21*t14*t23*t28-0.051897222222222*
        t1*t2*t3*t4*t8*t17);
    res[12] = 0.5*(-1.0567688679228883E-8*t1*t2*t27*t3*t31*
        t32/pow(rhoa,10.33333333333333)-0.001872433643323*t1*t13*t15*
        t19*t2*grada+2.366148652943867E-6*t1*t2*t21*t18*t23*t29);
    res[13] = 3.962883254710831E-9*t1*t14*t2*t3*t31*t32/pow(rhoa,
        9.333333333333334)-6.654793086404628E-7*t1*t2*t21*t5*t23*t30+
        2.340542054153986E-4*t1*t13*t15*t2*t20;

}

static void
revpbe_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[14];
 
    revpbe_fourth_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];

    ds->df3000 += factor*res[5];
    ds->df2010 += factor*res[6];
    ds->df1020 += factor*res[7];
    ds->df0030 += factor*res[8];

    ds->df4000 += factor*res[9];
    ds->df3010 += factor*res[10];
    ds->df2020 += factor*res[11];
    ds->df1030 += factor*res[12];
    ds->df0040 += factor*res[13];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        revpbe_fourth_helper(dp->rhob, dp->gradb, res);

    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

    ds->df0300 += factor*res[5];
    ds->df0201 += factor*res[6];
    ds->df0102 += factor*res[7];
    ds->df0003 += factor*res[8];

    ds->df0400 += factor*res[9];
    ds->df0301 += factor*res[10];
    ds->df0202 += factor*res[11];
    ds->df0103 += factor*res[12];
    ds->df0004 += factor*res[13];

}
